
import sys
import indexing.featureindex
from summaryprocessor import SummaryProcessor, SinglePositionProcessor, SkipZeroes
from genome.genomereader import GenomeReader

class FragmentCounts(SinglePositionProcessor):

  def __init__(self):
    SinglePositionProcessor.__init__(self)
    self.cts = {}
    self.totalFragments = 0
    self.fragmentPositions = {}
    self.genomicBase = None
    
  def getFragmentCount(self, base):
    return len(self.fragmentPositions.get(base, {}))

  def incr(self, base, seqpos, details):
    SinglePositionProcessor.incr(self, base, seqpos, details)
    key = base.upper()
    fragment = (details.spos, details.seq)
    self.cts[key] = self.cts.get(key, 0) + 1
    if not self.fragmentPositions.has_key(key):
      self.fragmentPositions[key] = {}
    if not self.fragmentPositions[key].has_key(fragment):
      self.totalFragments += 1
      self.fragmentPositions[key][fragment] = True

  def __str__(self):
    st = ""
    if not self.total: return ""
    bases = []
    for base in ['A', 'C', 'G', 'T', 'N']:
      if self.getFragmentCount(base)>0:
        bases.append((-self.getFragmentCount(base), base))
    bases.sort()
    for (negativeCt, base) in bases:
      if st: st += ","
      st += "%s=%s" % (base, -negativeCt)
    return st

def process(query, indexFilename, skipZeroes=(SkipZeroes.LEFT | SkipZeroes.RIGHT)):
  processor = SummaryProcessor(FragmentCounts)
  return processor.process(query, indexFilename, skipZeroes)

def printUsage(args, usage, minArgs=2):
  if len(args) < minArgs or args[1].startswith("-h") or args[1].startswith("--h"):
    print usage
    sys.exit(1)

usageString = """
Usage: python %s INDEX_FILE CHR MIN_POS MAX_POS
"""

chrs = ["chr"+str(x) for x in xrange(1,23)] + ["chrX", "chrY"]

if __name__=="__main__":
  printUsage(sys.argv, usageString % sys.argv[0], 3)
  if sys.argv[1]=="-ALL":
    queries = ((chr, '?', 0, 1000000000L) for chr in chrs)
    indexFilename = sys.argv[2]
  else:
    indexFilename = sys.argv[1]
    qchr = sys.argv[2]
    qmin = long(sys.argv[3])
    qmax = long(sys.argv[4])
    query = (qchr, '?', qmin, qmax)
    queries = [query]
  printedHeader = False
  for query in queries:
    printedInitialStep = False
    (qchr, qstrand, qmin, qmax) = query
    reader = GenomeReader(qchr)
    prevPos = 0
    for (x, cts) in process(query, indexFilename, SkipZeroes.ALL):
      genomic = reader[x]
      if x != prevPos+1 or not printedInitialStep:
        print "fixedStep chrom=%s start=%ld step=1" % (qchr, x)
        printedInitialStep = True
      prevPos = x
      if cts.getFragmentCount(genomic) != cts.totalFragments - cts.getFragmentCount('N'):
        # print "%s\t%ld\t%s\t%ld\t%s" % (qchr, x, genomic, cts.total, str(cts))
        print "%ld\t%s\t%s" % (cts.totalFragments, genomic, str(cts))
      else:
        # print "%s\t%ld\t%s\t%ld" % (qchr, x, genomic, cts.total)
        print "%ld\t%s" % (cts.totalFragments, genomic)

